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Abstract 

This paper investigates the solution of a paraboHc inverse problem based upon 
the convection-difFusion-reaction equation, which can be used to estimate both 
water and air pollution. We will consider both known and unknown source 
location: while in the first case the problem is solved using a projected damped 
Gauss-Newton, in the second one it is ill-posed and an adaptive parametrization 
with time localization will be adopted to regularize it. To solve the optimization 
loop a model reduction technique (Proper Orthogonal Decomposition) is used. 

Keywords: Inverse problem, regularization, adaptive parametrization, time 
localization. Finite Element method Proper Orthogonal Decomposition 



1. Introduction 

Inverse heat or mass convection problems, classically deal with the estima- 
tion of wall heat flux densities or intensities of source terms [321 HZl ttH HI] ■ As 
mentioned in 11], inverse problems are usually mathematically ill-posed and 
regularization methods have been developed to ensure stable solutions. For an 
overview, see [33[ I15[ I17j for example. Classical methods are penalization such 
as Tikhonov's regularization |32| . or Bayesian methods using prior information 
f7 , iterative regularization and regularization using singular value decompo- 
sition followed by truncation of the singular values spectrum [| BD] . 

In this paper we are interested in solving an inverse convection problem, 
whose direct model coincides with a parabolic convection diffusion reaction equa- 
tion on a fixed domain. To deal with its ill-posedness we adopt a regularization 
algorithm based upon Truncated Singular Value Decomposition (TSVD) and 
diagonal scaling [26 ; moreover an adaptive parametrization with time local- 
ization is formulated, to reduce the computational cost of the Gauss Newton 
algorithm and to obtain a better conditioned sensitivity matrix (cfr. e.g. figure 
14]). 

Convection-diffusion-reaction equation can be used to model a variety of 
physical problems. For example in [9 , this equation is used to predict water 
quality in rivers, by measuring the quantity of organic matter contained. The 



importance of these pollution is estimated by the measures of the so-called BOD 
(Biologic Oxygen Demand) and COD (Chemical Oxygen Demand). In [3] the 
problem of identifying the location and the magnitude (intensity) of pollution 
point sources from the measurements of BOD on a part of the river is con- 
sidered: the problem of source term identifycation is solved using an algorithm 
based on the minimization of a cost function of Kohn and Vogelius type. Also in 
|29j water pollution is considered: knowing the origin of the source of contami- 
nation is probably the most important aspect when attempting to understand, 
and therefore to control, the pollution transport process. Thus, a challenging 
issue in environmental problems is the identification of sources of pollution in 
waters. ^29] deals with source identification problem, using Boundary Elemet 
Method (BEM). In [TT], the same problem of source estimation is considered 
to estimate the time-varying emission rates of pollutant sources in a ventilated 
enclosure, assuming that the velocity field is stationary: in fact in the frame 
of occupational risk prevention, the knowledge of both space and time distri- 
butions of contaminant concentration is a crucial issue to evaluate the workers 
exposure. Althought air pollution is considered, instead of water, the under- 
lying model is still a convection-diffusion-reaction equation, with a different 
convective velocity field. In [11] source's location is supposed to be known. Pos- 
sible applications of this study are concerned with cartography of pollutants in 
buildings, estimation of contaminant emission rates inside manufactures, leak 
detection, environment and process control through 'intelligent sensors' (con- 
trolled ventilation with closed-loop function of pollution threshold). A similar 
problem is considered in [5]. Finally in [7], a convection inverse problem is 
solved to determine an estimate of the source term as a funcition of the altitude 
and the temporal of iodine-131, caesium-134 and caesium-137 in the Chernobyl 
disaster. 

In general, in inverse convection problems, either distributed control |llll29j . 
or boundary control [34J or both [20] are considered. In the present paper we 
are interested in estimated location and intensity of pollution, and we assume 
to deal with boundary control, i.e. we suppose that the sources are located 
along domain's boundary. Thus, as in |34j . we deal with an inverse problem in 
which one is looking for the unknown conditions in part of the boundary, while 
overspecified boundary conditions are supplied in another part of the boundary 
(here the outflow region). As mentioned above, this type of problem can model 
both water and air pollution. 

As mentioned e.g. in |13| . in inverse problems or optimal control or opti- 
mization settings, one is faced with the need to do multiple state solves during 
an iterative process that determines the optimal solution. If one approximates 
the state in the reduced, fc-dimensional space and if k is small, then the cost 
of each iteration of the optimizer would be very small relative to that using 
full, high-fidelity state approximations. Thus Proper Orthogonal Decomposition 
(POD) will be adopted in this paper as model reduction technique, to bring our 
study closer to a real time problem. 

In section[2]the direct problem is described. In section[3| the inverse problem 
is formulated. Section [5] deals with the problem of known source location, while 
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in section [6] also source position is estimated. 

2. Description of the direct problem 

Let [0, tf) gM. and O be an open, limited and Lipschitz continuous boundary 
subset C M^, sufficiently regular. We denote with dH. the boundary of ft. Let 
C : [0,tf) X fl ^ R, C — C{t,:x.) be the solution of the following (direct) 
parabolic convection-diffusion-reaction equation: 

in {0,t-f) X fl 
on {0} X fl 

on (0,t/) X Tin (1) 
on (0,i/) X r„p 
on (0,t/) X Tdown 
on {0,tf) X Tr 

where Ti„, r„p, Tdown and are given disjoint sets such that dfl = TinUTupU 

^down U . 

Suppose that dn e H'^^Tin), Cup € H^iTup), the initial condition Cq € 
L^(f2) and the coefficients are independent on time, moreover G 

> Mo > for all X e f], cr G a{x) > a.e. in fJ, u G [L°°{Vl)Y, 

div(u) G L^(f2) are known. The direct problem consists in finding the concen- 
tration C over f2 at time tf. 

As in , we assume that the physical properties of the fluid are constant 
and that the transported contaminant is considered as a passive scalar, which 
means that it does not affect the velocity field. Thus we suppose to know u. 

An example of the 2D domain i7 is illustrated in figure [T] 
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Figure 1: Example of problem's domain Q. 

2.1. Wellposedness of the direct problem and finite element discretization 

Let iJp ur (^) t>e the set of G 77^ (il) s.t. t;|^ = 0. Given 

.ur„„uri„(^)' weak formulation of (|l|) consists in finding C G 
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(2) 



L^{o,tf;HHn))nc''i[o,tfy,L^{n))s.t. 

iiC{t),v) + a{u{t);C{t),v) = 0, e V, 

C(0) Co, in n, 

where a(u; •, •) is a bilinear form defined as 

a{u;w,v) := / kWwWvdu + / u ■ Wwvdu + / awvduj. 
Jn Jn Jn 

The wellposedness of the variational formulation is studied e.g. in [28]. 

Consider now two families of subspaces {Wh, h > 0} and {Vh, h > 0} of 
H^{il.) and V respectively, and let Cq^/i S be a suitable approximation 
of Co. Then the Finite Element (FE) discretization of ([3| consists in finding 

Ch e Wh s.t. 

Ch{0) - Co,,., inn. 

Given a basis of {(/i,;}! * = 1, . . . ,Nh, where Nh denotes the number 
of nodes in fi, it is well known that the FE discretization is equivalent to the 
solution of the following system of ODE's: 

MCit) + A{u{t))C{t) - F(C„), 

C(0) - Co. ^ ' 

where Mij = ^(u)^^ = a{u; (j>i, (f>j) and r(Ci„) involves boundary con- 

ditions, in particular Ci„. 

Given a time step Ai, consider a uniform subdivision of [0,^/) s.t. (TV — 
l)At = tf. Discretizing Q in time, using e.g. the implicit euler method, we 
obtain 



{M + AtA{u{k+l)))C{k + l) = MC{k) + AtFiCn), 

C(0) - Co. 



(5) 



2.2. Proper Orthogonal Decomposition (POD) reduction 

To obtain a faster solution algorithm, we adopt a reduction technique. A 
complete overview of all classical methods can be found e.g. in [HEI]. Since the 
right hand side in Q depends on boundary conditions, it varies at each iteration 
of the optimization loop used to solve the inverse problem. As a consequence 
techniques largely used for linear constant matrices problems, like e.g. Balanced 
Truncation (BT), becomes too costly to be used. Thus we choose to adopt the 
Proper Orthogonal Decomposition (POD) method: althought its basis is stricly 
related to local dynamics, it is less costly to compute. In this paper we are 
focusing on the inverse problem solution strategy, thus we will not enter in 
details in the description of POD, we only summarize the main aspects: the 
interested reader can found a complete overview for example in [181 114] . 
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Given a time step At > (which could be different from At), consider 
tm G (0,t/) and N s.t. TV At = tm- first the unreduced model (Wl) is solved 
in [0,t„i], collecting the matrix of snapshots X — (Cj), where Cj G K^'' is 
the nodal vector of the FE discretization at time tj = JAt, j = 0, . . . ,N. 
After computing the Singular Value Decomposition (SVD) of X, X = USV* , a 
suitable threshold k is chosen. A largely used strategy is to choose k s.t. 



is greater than a fixed toUerance. Another possibility is to impose that the first 
k singular values are greater than a fixed tollerance > 0. 

It can be proved [M] that the k-th POD basis {ui}, i = 1, . . . , k, Ui U{: 
, i) € M^'' is the solution of the following minimization problem 



N 

min y 



k 2 

1=1 



(6) 



i.e. for every fixed k the mean square error between the elements C{tj) and the 
corresponding k — th partial sum of ^i=i{C{tj) ■ is minimized on average. 

Finally Q is projected on the space generated by the first k POD basis 
vectors, i.e. we solve the reduced system 



UlMUka{t) + UlA{u)UkaL{t) = C/*F(a, 
a(0) = C/*Co. 



(7) 



in {tjn,tf), where Uk ■= U{:,1 : k), i.e. the system is projected on the subspace 
generated by the first k columns of U. We denote with C{t) := UkSi{t) the 
estimate of C{t), t G {tm,tf) computed using POD. 

3. Continuous inverse problem formulation 

We are interested in solving the following inverse problem: given the addi- 
tional a priori information 

C = Cs, on [0,tf) xVdown, (8) 

where Cg € C'^{{0,tf), L'^{Tdown)) is a known scalar function, determine C,*„ € 
Hi(rin) such that 

Qn = arg mill -Fd(C„0. (9) 

c,„e//2 (ri„) 

where the cost function is 

MC^n) ||C(a„;i,x)-C,(t,x)||^.([o^,^]^p^_) =^ {C{an;t,^yCsit,^)fdjdt 
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and we have explicited the dependence of C, sohition of Q, on Cm- C(Ci„; t, x) := 
C(t,x) s.t. C(t,x) = Ci„{x), when x G . 

As mentioned in [34] , one may consider Cg to be a desired one. In that case, 
the present inverse problem is a design problem where the boundary flux Ci„ is 
controlled such that a desired concentration is achieved on the boundary Tiiown- 
Cs can also be considered to represent a continuous approximation of a set 
of discrete experimental temperature measurements obtained at finite number 
of locations in the boundary Tdown and at discrete time instances within the 
interval [0,^/). In this paper we mainly refer to the second case. Observe 
that this class of inverse problems are of significant experimental interest for 
situations where the direct measurement of the heat flux Cin is not possible. 

As indicated in |34) . the main difficulty with the minimization problem (|9| 
is the calculation of the gradient of J-. Mainly two different approaches could 
be used: the first discretize than optimize or vice versa the first optimize than 
discretize. In this paper we focus on the first strategy: in particular we adopt a 
discrete approximation of J- (Ci„), combined with a Gauss- Newton approach, 
as explained starting from section [5] 



4. Formulation of the discrete inverse problem 

In the first discretize than optimize context, Cs(i,x) is known only in the 
Uy nodes of Tdown, for every discrete time tj, j = 0, . . . , N — 1. We denote with 
Cs{j) G K"" the vector of measured concentration at iteration j. 

As in [llj , we assume that the flow dynamic boundary conditions are steady 
state and for simplicity we suppose that 

ng 

r - 1 Ir^'^ 

1=1 

being F-^ disjoint sets, such that Cm is constant on each T^^, for all I = 
l,...,ne. 

Thus we have to estimate a vector i9 of ng non negative parameters: equiv- 
alently we assume that the function Cm G H^iTm) could be identified by a 
piecewise constant function Cm{'&) such that 

Cmm^)^m. xe T^i 

Since we are solving an inverse problem we indicate with ■d the estimate of the 
realnarameters Let 11 : M^'* M"" be the map which projects the solution 
of (k| on the Uy nodes of Tdown'- thus we denote with n(C(i^;t)) the estimate 
of C(Ci„(i?); x) at time t on Tdown {predicted concentration) obtained solving 
^ imposing Ci„(i?) on Tm- 

In a space-time discrete setting (|9| could be restated as 

i9* = arg min Fd{^), (10) 
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where the discrete cost function is defined as 



N 



-^•^W '-^nT. l|n(C(Q„(,9); j)) - CsiM: 



(11) 



4.1. Proper Orthogonal Decomposition reduction 

Using model order reduction techniques to solve ([9]), consists in replacing 
the cost function (111 in (10) with the following one 



1 ^ II 



(12) 



where C is the solution of ([7|). An example of application of POD to solve 
optimal control problem can be found e.g. in |14) . 

Since the POD basis depends on the collected snapshots, it is necessary to 
update the projection space as the estimated control Cm varies. Let n a small 
positive integer: at every iteration i in this paper we adopt the following index 



^C'(a„(#'));j)-C(C„(#'));j) 



i.e. we compare the first iterations of the unreduced system with those obtained 
projecting on the old POD basis used at iteration i — 1. Only if I^*) is greater 
than a fixed threshold, the i-th basis is updated, computing new snapshots, as 



described in section 2.2 Two strategies can be used [Hj: old snapshots can be 
discarded or not. In practice this consists in adding POD modes computed in 
the i — 1-th iteration to the new snapshots ensamble: in this case the projection 
space is more robust to control variations but usually is slightly bigger. For our 
experimental tests we prefer to discard old snapshots. We obserse that in [l4] 
a new basis is computed at every iteration, without considering an index I. 



5. Known source location Tin 

As a first step toward the solution strategy, we consider a simpler problem, 
assuming that the source location Tin is known. 

5.1. Solution uniqueness 

In this section we demonstrate that if Tin is known, then the discrete inverse 
problem admits a unique solution, since there are no local minima. Moreover 
changes in Ci„ corresponds to changes in the registered concentration. 

First of all we prove the following Lemma, which justifies mathematically 
the physical principle that, as Ci„ increases on Tin, the concentration on Tdown 
increases too. 
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Lemma 5.1. Consider the two problems 
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represented in Figure ^ (up ), where 
for every x G Tin. 

Then C2{t,x) > Ci{t,x) for every t G (0,t/) and xE r^, 



(13) 
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Figure 2: 



Proof. 

Define w := C2 — Ci, which solves 



dw 
dt 



fiAw + V 



(uw) + (7W 
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= 0, 


on 


{0} X n 


W 


_ ^(2) ^(1) 


on 


{o,tf) X r„, 


W 


= 0, 


on 


{o,tf) X r„p 




= 0, 


on 


(0, tf) X Tdown 


w 


= 0, 


on 


(o,t/) X r. 



(14) 

as illustrated in figure [2] (down) . 

Observe that w is smooth only inside the domain, but it is not continuous 
near the boundary, where it admits discontinuities of the first kind: thus general- 
ized solutions must be considered. The strong minimum principle for parabolic 
operators can be extended for generahzed solutions [TUlIin] (cfr. appendix ??): 
thus we know that the minimum is assumed at the boundary. 

Moreover, for every open neighbourhood U of Tdown, such that w is regular 
inside U Hfl, ^ = on Tdown implies that the maximum and the minimum of 
w over U nil cannot belong to Tdown (cfr. [TO]). 
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As a consequence w > in (0, t/) x and, since the minimum is not attained 
on Tdowm w = C2 — Ci > on Tdown, for all t £ (0,t/) i.e. the thesis holds 
true. 

□ 

The following Proposition is equivalent to prove that there are no local min- 
ima. 

Proposition 5.1. For every '& E M"", i9 7^ there exists at least a sequence 
of profiles — '^^ converging in C'^iM."'") to the real profile "d* , such that 

Proof. We can construct the sequence {■»?«}„ in the following way. For 
every k = 1 . . . . , 

Thus 'dng = 1?* by construction. Moreover the corresponding sequence of cost 
functions is decreasing: Tdi'&i) > J^d{'&2) > ■ ■ ■ > J^dii^*)- This fact is a direct 
consequence of the application of Lemma [5.1[ suppose that -dk-iik) < •d*{k). 
Then 'dk{k) > dk-i{k) by construction and thus n(C(i5fe; t)) will be higher than 
n(C(i9fc-i;i)) for every t € (0,i/) (Lemma [slj) and thus closer to n(C(i9*;t)). 
Analogously if -dk-iik) > '&*{k), applying Lemma 5.1 n(C(i?fc; t)) will be lower 
than n(C(i9fc-i;i)) for aU t and thus closer to n(C(i?*;t)). 

□ 

5.2. Numerical solution strategy 

Starting from an initial guess \ line search algorithms find the fc + 1- 
iteration starting from the fc-th one in the following way: 

where the damping parameter a^'^) is obtained using a bisection procedure. 
The standard Newton step consists in solving at each iteration the system 

Let TZ : M"*^^ W^y^ be the map s.t. starting from an x matrix 

/ h 

B = [61, ... , 6jv], it gives 7^(5) := 

V hN 

The computation of J-^ , Hessian of the cost function, usually is expensive. 
Moreover, since we are dealing with a least squares problem ([9]), we adopt the 
Gauss-Newton approximation (cfr. [H]), i.e. we solve 

vf^^.jsW =e^(.,, (16) 
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where the sensitivity matrix ^^c^) G jg such that 

^.c^oO,*) -:|^7^(^(C(#V))), 

for all i — 1, . . . , ne and the prediction error 

e^,., :=7^(C,(.))-7^(^(C(#V))). 

To compute numericahy the sensitivity matrix a finite difi'erence scheme is 
needed: 

V j) « ^ [7^(^(C(,§^'^ (1), . . . , S^"^ (j) + J, . . . , 1?^'^ K); •))) - 7^(^(c(19^'^ ;•)))" 

where 5 > is a smaU perturbation parameter. 

Observe that in general this approximation is computationahy expensive, 
since, it requires the computation of the concentration also for the perturbed 
input. When is known, only very few parameters are considered, thus 
this approximation is effective. The problem becomes more involving when Tin 
is unknown, since the number of parameters is higher: in section |6] we will 
explain how the adaptive parametrization and time localization can reduce the 
computational cost. 

If (5 > is too small, the finite difference estimate could be inaccurate, since 
at the numerator we are considering the difference between two quantities which 
has approximately the same absolute value, and this is divided by a very small 
denominator, which amplifies the error. A possible solution e.g. is to adopt the 
Complex-Step Derivative Approximation [23] . in which an immaginary incre- 
ment iS is used, approximating 

*^(.)(:, j) « ]lm (7^(^(C(1?*'\l), . . . , +iS,..., ^^'\ne); •)))) . 

Finally observe that we are assuming that the pollutant is put into the 
domain, thus 

> : 

as a consequence we need also a projection step onto [0, +oo) of each component 
~ (fc) 

of 1^ , after its computation. 

5.3. Numerical results 

In this section the Projected damped Gauss Newton is compared to other 
classical solution strategies. Experimental data are simulated numerically, on 
n = [0,8] X [0, 1], Th = [0,8] X {1} U [0,8] x {0}. Moreover the velocity field u 
is modelled as a Poiseuille flow i.e. 



u(a::i,X2) 



— AVX2 + 4:1^X2 





10 



We assume that v = 50, ii — 0.1, a = 0.1 and Cup = 0.1. Moreover in this 
section a Gaussian error of variance 0.05 and mean zero is added. 

Classical solution strategies cited in this section are well described e.g. in 
|16j . As a regularization parameter, when needed, we use a = 0.01, moreover 
we choose a maximum number of iterations maxn = 20. Consider the following 
two sparse examples: 

1. r„, = [4,4.5]x{l}, i? = 100; 

2. r,„ = [4.5, 5] X {1} U [1.5, 2] X {0} , i? = (100, 80); 



and see how different techniques approximate them. 

First of all we consider the example 1. Performances of different methods 
are depicted in figure [3j In the second example, two parameters have to be 
estimated: results are plotted in figure |4j 

Observe that in both cases the Projected damped Gauss Newton algorithm 
performs well, converging faster to the optimal solution. It should be noted 
that, in constrast to Tikhonov and Levenberg-Marquardt it does not need a 
regularization parameter: this is important because it tells us that, knowing the 
source location, the inverse problem is not ill-posed, as stated in Proposition 



5.4- Reduce the order of the system using POD 

In this section we analyze the POD reduction introduced in section |2.2| on 
a test case. Consider example 2 introduced in the previous section; in POD re- 
duction two parameters plays a central role: t™, which characterizes the interval 
[to, trn] when snapshots are collected, and the threshold on the singular values 
of the snapshots matrix. As can be seen in table [T| increasing i„ corresponds 
to a better approximation, since more snapshots are collected. To obtain higher 
accuracy decreasing tm, it is necessary to increase To-, i.e. to consider a bigger 
reduced model. 
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Table 1: Example 2 of section 5.3 choosing different intervals [to,im] to collect snap- 



shots and different thresholds To- on singular values of the snapshots matrix. 
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Figure 3: First example: different strategies. Left: cost function and error, right: 
convergence. First row: projected damped Gauss Newton, second row: Levenberg Mar- 
quardt, third row: steepest descent, fourth row: Tikhonov method. 

It is important to note that the reduction is significative with respect to the 
unreduced model, which has dimension 1071. However, as described in section 
|2.2[ it should be noted that it is necessary to update the POD basis: in all 
these examples the basis is updated at every new iteration, imposing 0.1 as a 
threshold on I*^*^ . 

A more involving problem is considered in section |6j where it is assumed 
that also the source location Tin is unknown. In general in that case pro- 
jected damped Gauss Newton could not be sufficient and it is too costly, thus 
it is necessary to adopt a suitable solution strategy based upon an adaptive 
parametrization and time localization. 



12 



Figure 4; Second example: different strategies. Left: cost function and error, right: 
convergence. First row: projected damped Gauss Newton, second row: Levenberg Mar- 
quardt, third row: steepest descent, fourth row: Tikhonov method. 

6. Unknown source location Tin 

Suppose now that the location of Tin is unknown. 

6.1. Ill-posedness of the problem 

To study analytically what happens when Tin is unknown, we consider a 
simplified model problem: let C = C{x), x € [a:i,a:2] C M., X2 > xi be the 
solution of the following one dimensional ODE: 

— /iC (x) + uC (x) ~ f{x)^ in {xi,X2), 

C{xi) = Cup, (17) 

C'{X2) = 0, 

-h-^ - { 1; il;Jrele ix,,X2) ' ^ > 0' e (---^)' ^ ^ 1) 
s.t. x„i zt h E {xi,X2). Observe that (|17[) can be viewed as the one dimensional 
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stationary counterpart of ([T]) when a = and considering only the x-axis in fig- 
ure [ij the unknown immision boundary Tin can be represented by an unknown 
forcing term /, apphed in [xm — h, Xm + /i], of intensity M. In this context the 
inverse problem (|9| is equivalent to determine the source position {h and x^) 
and intensity (A/) given the measured concentration Cs G M in X2. 



Problem (17) can be solved analytically, obtaining 



Ci+C2ei^^, x<Xm — h, 

C(a;) = <( C3 + ^x + C4e>, \x - x,n\ < h, 
cs + ceCM^, x>Xjn + h, 

where Ci , . . . , cg are suitable real coefficients obtained imposing boundary con- 
ditions and continuity of u and u in Xm ± h. In particular we are interested in 
estimating the concentration at the measurement point x — X2- For simplicity 
we assume that xi — and X2 — 1- Thus it can be derived that 

„ 1 -u{x„i + h)f u{x„i + h) ( 2uh 
cg = 0, C5 = ^ exp zunM exp h /iM 1 — exp 

Thus C{x) is constantly equal to C5 in [xm + h,l]. We can now study how 



Figure 5: Solution of {17) at the measurement point X2 ~ 1 for different values of M 



(left), h and Xm (center), Xm (right). 

C(l) depends on M, /i, L. We consider fi — 0.5 and u — 10 (Peclet number 
Pe = 5^ = 10, quantity that characterize convection diffusion problems). As 
can be seen in figure [5] varying only M, fixing h and Xm (i-c knowing the 
source location), corresponds to a linear striclty increasing C(l) (cfr. figure |5] 
(left)). On the contrary fixing M but varying h and Xm corresponds to the 
surface plotted in figure [s] (center) : fixing h for different values of Xm we obtain 
almost the same C(l) (cfr. figurejs] (right)). Thus measuring C(l), the problem 
of determining the source is ill-posed in the stationary regime. Increasing the 
Peclet number this phenomenum is stressed. 

Even for this simplified ID stationary problem, in general unknown source 
position gives rise to an ill-conditioned problem. 

6.2. Numerical solution of the discrete inverse problem 

The problem consists in estimating both the position of the sources Cm in 
the horizontal segments Th := Tr U Tin and their intensity. 
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6.2.1. Algorithm 1: working on the finest subdivision 

First of all we consider < xi, . . . , a; (/> > a reference uniform finest sub- 



I Si Ul I 

division of Th of step length Ax, which represents the minimum amplitude of 
estimated source emissions. 

The simplest strategy consists in applying the Gauss Newton method di- 
rectly on the finest subdivision (cfr. algorithm [l]) , i.e. in estimating n^f^ pa- 
rameters. This problem is particularly demanding for its high computational 
cost, due to the large number of parameters to be estimated at each Newton's 
iteration. Moreover it should be noted that if we want to estimate a sparse 
vector of parameters, working only on the finest subdivision is not efficient, as 



we will see in the following sections. To solve the system (16 1 both TSVD and 



Algorithm 1 Sketch of the algorithm working on the finest subdivision: 

^0 



Given the finest subdivision of F^, = 0, f/* = 1; 
while ) < tol do 

solve ipgks'^ = egk ; 

= e'' + ^''s" 

projection: for every j G [0,'ng — 1] s.t. 9''^^(j) < 0, impose 0'''^^{j) 
compute 

if > then 

I = 0; 

2 

while Tn^{e^*^) < Fn{0^) do 



11: 0*'+' =0'''+m'=''s''' 

12: 1 = 1 + 1; 

13: m"-' = ^ 

14: end while 

15: end if 

16: end while 



diagonal scaling [26 arc used. The last one, presented in [21] to solve a con- 
duction inverse problem, works as follows: at iteration k, given the subdivision 

iSC^) = \ xi, . . . ,x (/) k for every i = 1, . . . , n^c \ where ni'^^ = n^P denotes 



the number of columns of ^^C') iteration fc, ^^(fe)(:,j) is multiplied by a 

weight dj, equal to the length of the maximal segment of the current subdivi- 
sion, divided by the length of the segment corresponding to the i-th column. 
Thus diagonal scaling corresponds to solve 



max 



vP^c.ii^WiW = e^(.,, D^^) =diag{df^), - ^2±l^ 



(18) 



instead of ( 16 1. 
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6.2.2. Algorithm 2: working on the finest subdivision with time localization 

As explained in section [6^] in the stationary regime the problem is illposed: 
time localization corresponds to a better conditioned problem, since it consists 
in selecting only those rows of the sensitivity matrix which are significative for 
the dynamic, i.e. corresponding to the transitional dynamics. 



Figure 6: Example of partition of Q in sections. 

More precisely, the idea is to partition the domain 17 in a suitable number 
> 1 of sections lA = {sj}, j = 1, . . . , (cfr. e.g. figure[6]). Referring to figure 
[l} we suppose that Sj := [ij^ij+i] x [2/1,2/2], ^1 = xi, Cn.+i = X2. In particular 

in algorithm 2 we assume that {^i,...,^„ +1} = <a;i,...,a; (/) >, i.e. it 

coincides with the finest subdivision. Denote with the set of parameters 
belonging to Sj . 

The algorithm works as follows: starting from s„^ , it computes the sensitivity 
matrix only of those parameters belonging to only in the time interval 

[4"°'*'^/"°^]' 4"°'' — ^0, ^ if- below it is explained how to choose the 

interval. The estimate of the parameters of section is done as explained 
before, using a projected damped Gauss-Newton method. To regularize the 
problem both TSVD and diagonal scaling are used. 

Define O^-'-' , j — 1, . . . , — 1, the set of parameters estimated in section Sj+i 
greater than a threshold €3 > 0: then in section Sj all parameters belonging to 
qU) u jU) -^vill be estimated, only in the time interval [ijj''', In algorithmji] 
previous ideas are summarized. 

6.2.3. Algorithm 3: using an adaptive parametrization 

The alternative is to use an adaptive parametrization, i.e. to adaptively 
update the subdivision of Th used in the current iteration of the Newton method. 
This strategy is particularly indicated when dealing with a sparse vector of 
parameters: in this situation it is important to localize T^^ in Tft and to refine the 
parametrization possibly only around that point. It reduces the computational 
cost reducing the number of columns of the sensitivity matrix. A similar strategy 
has been presented in [S], to solve an inverse conduction problem of corrosion 
estimation. 

The algorithm works as follows: starting from an initial coarse subdivision 
of r^, iS^^-*, at the k-ih. iteration the algorithm first computes a Gauss-Newton 
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Algorithm 2 Sketch of the algorithm working on the finest subdivision with time 

localization: 

1; Given {gi, • • ■ , Cris+i} coincident with the finest subdivision of and the threshold es > 
0; 

2: while Tdi'^') < tol do 
3; i = ns; 0^"=' =0 
4; while i > do 

5: Let be the set of parameters of S that belongs to section i; 

6: in [to^'it^''] apply the regularized projected damped Gauss Newton method to opti- 

mize parameters whose indices belong to U O''^; 
7: update the positions U O**' of ; 

8: define 0^'~^^ as the set of indices of parameters greater than £3; 

9: i = i - 1; 

10: end while 
11: end while 



iteration -0 G M"" . For every element of 1? G M"" greater than a fix thresh- 
old ei > 0, the segment of 5^'^'^ corresponding to that parameter is bisected: 
thus a new subdivision S'^'^^^^ is defined adding to 5^*''-' all the computed mid- 
dle points. Finally are selected only those parameters which are greater than 
a fixed threshold £2 > 0, and we indicate with A*^'') this ensemble; the other 
parameters remain constant in the following iteration. The main ideas of the 
adaptive algorithm are sketched in algorithm [3j 

To avoid a large over-refinement, the bisection procedure can be limited, for 
example applying it at each iteration only a certain number of times, choosing 
the segments to be refined as those corresponding to greater parameters. 

6.2.4- Algorithm 4-' using an adaptive parametrization and time localization 

The idea now is to reduce both the number of columns and of rows of the 
sensitivity matrix. As in algorithm 2, the domain is partitioned in ng > 1 
sections U = {sj}, j = l,...,ns, however in algorithm 4 we assume that 
{Ci: ■ • ■ J Cris+i} = iS'^-*, i.e. it coincides with the coarse initial subdivision 
applied in the adaptive strategy. In section s^, considering the time interval 
[^q"''', <j "*], all parameters belonging to O'-'^ U 1^^^^ will be estimated, and the 
adaptive procedure will be applied until a minimum is reached. Observe that 
this coincides with an internal loop: the ideas are summarized in algorithm |4] 

6.2.5. Time localization: how to choose time intervals [t^Kt^^^] 

A key point is the choice of the local time intervals [tp''*, i^''], for every 

(i) (i) 

section Si, i = l,...,ns, > to and ty < tf. The i-th interval must be 
chosen such that it contains the transitional dynamics of section Si but not 
that of sections Sj, j < i. To describe more clearly this idea, consider the 



model problem introduced in section 5.3 moreover suppose for simplicity that 
Ug = 2, {^1,^2,^3} = {0,4,8}, as depicted in figure [7j and consider as the 
finest subdivision a uniform one of step length 0.5. Consider figure [Hj the j-th 
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Algorithm 3 Sketch of the adaptive algorithm: 



1: Given the finest subdivision of r^, of step length Aa;, the tolerance tol > and thresholds 
ei, £2 > 0, consider the coarse subdivision S^^^ = < x\, . . . i >, of Fh.; 

2: 1?^ = 0„i e K"9; 

3: / = 1, A(i) = [1, . . . , rig], set of indexes of parameters to be optimized 
4: while J^d(^') < tol do 



6: 
7: 
8: 
9: 
10 

11: 
12: 
13 
14: 

15 
16 
17: 
18 
19 
20 
21 
22 



for all i e [1, /] do 

if 0^(i) > ei% bisect the corresponding segment then 

n^+^ =n^+^ + 1, / = J + l; 

let [a;'+^(6'(i)),3;'+^(9'(i))] be the segment corresponding to parameter 
50+1) ^ u 
end if 
end for 

given the subdivision 5''+^^ apply the projected damped Gauss-Newton method, op- 

timizing only parameters whose indexes belong to A*^^ , obtaining i? S K"e 

AC') = 0; 

for all i e [1, /] do 

if > £2 then 

A(*) = At'') Ui; 

end if 
end for 

1 = 1 + 1; 
end while 
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Algorithm 4 Sketch of the adaptive algorithm with time localization: 
1: Given the partition of H {^i, . • • ,^ns 

+1} = 5(1), the thresholds ei, €2, es > 0, = 0, 

i = 0; 

2: while Td{'^^) < tol do 
3: j=j + 1; 
4: i = n^; ©("^'l' = 
5; while i > do 

6: ' = li A^'"^^ = [1, . . . , raj''^], set of indexes of parameters to be optimized 

7: while a minimum is reached do , ^ 

8: Let be the set of parameters of that belongs to section i; 

9: in [tg^^tj^] apply the regularized projected damped Gauss Newton method to 

optimize parameters whose indices belong to U 

10: update the positions U O^^''^ of ; 

11: apply the adaptive strategy: 

12: 50,i,(+l) = 5O'.».0| 

13: nj*''+' = n^'*^ 7 = ''+^ 

14: for all ifc 6 [1,7] do 

15: if eJ''''(fc) > ei then 

16: update 5^^''''+^', bisecting the segment corresponding to ^■'''■'(A;); 

17: end if 

18: end for 

19: given the subdivision S^-'"*''"'"^-' apply the projected damped Gauss-Newton 

method, optimizing only parameters whose indexes belong to aO'*''+^) 

20: \U,i,l+i) = 0; 

21: for all k e [1,/] do 

22: if t?J''''+l(A;) > €2 then 

23: aO'*''+i) = aO''>'+i) U k; 

24: end if 

25: end for 

26: i = i + 1; 

27: if the subdivision has been refined, update 

28: end while 

29: 5^^'''' = 

30: define (pCj''"!) as the set of indices of parameters greater than €3; 

31: i = i — 1; 

32: end while 

33: # = 1?-''*' 

34: end while 
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Figure 7: Partition of Q, m 2 sections to draw the curves of figure to obtain the 
red (blue) curve of figure^ it is considered the mean concentration on Tdown, ob- 
tained imposing a control different from zero only m the most left segment of the finest 
subdivision of the upper horizontal segment of section S2 (si), indicated in red (blue). 



curve Q, j = 1,2, represents the mean concentration (left) and its derivative 
(right) at the outflow when the boundary control is different from zero only in 
the most left position of Sj with respect to the finest subdivision. The interval 
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Figure 8: Time evolution of the mean concentrations at Taown, Ci '^^d. ^2, (isft) cind 
their derivative (right) for different boundary controls: the boundary control is different 
from zero only in the most left position of the finest subdivision of si (blue) and S2 
(red). 



corresponding to S2 can be [t};',t) '] = [180,260], when the red dotted curve 
corresponding to S2, C2 i is increasing (transitional regime) and the blue curve 
corresponding to si, (^1, is flat, i.e. when only the pollutant put into SI in 
could reach r^oii^n- While in si the choice can be = [240,400], since in 

this interval the transitional regime of si occurs, as showed by Ci. This intervals 
are used in section 6.4 to test algorithm 4. 

The previous idea can be extended more rigorously to a general number of 
sections: let Q, i — 1, . . . , n^, be the mean concentration at the outflow r^o^,„ 
when the boundary control is different from zero only in the most left position 
of Si, with respect to the finest subdivision. Consider a small threshold £4 > 0, 
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and two positive parameters d, D > 0. Given 



4"°^ = mintg[t„^f^^] |C,V,(0 > e4 and Cl^_i(t) < , 
[Cjt)>ei andC,_i(i) <e4}, 

then for i = 1, . . . , ti^ — 1 



maxtg[4^^t^] |Ci(i) > £4 and C--i(t) < £4} , « > 1 

min|i}+i + £>, maxigft^^^^.] |C-(t) > £4 and C'-iW < £4}} , i = 1 



The parameter d allows a small overlapping between local time intervals, 
while D could limit the length of the inteval [4^\^/^'']- ^.g. in the example 
presented above, considering — 2, d = 0.2 and D = 1.6. 

Observe that the definition of the intervals [io\^/^] depends on the shape 
of the domain, on the velocity field and on the coefficients of the PDE ([T]): each 
time one of them is changed, also the intervals should be estimated, observing 
the transitional dynamics of each section, as explained above. 



6.3. Comparing computational costs 

In this section we compare the computational costs of the four algorithms. 

The first one consists in using the finest subdivision, with the projected 
damped Gauss Newton strategy. The computational cost of each iteration is 
pretty high: since the solution of the direct problem has cost NN^, comput- 
ing the sensitivity matrix e jj"y^x"6i^' has cost n^g^NNf^, where n^^^ is 
the number of parameter of the finest subdivision, which is maximal. More- 
over computing the SVD to obtain the new iteration has cost An^N'^n^P + 
9)Nny{n^J^Y -\- 9{n'^/^Y . Finally computing the new prediction error has cost 

To decrease the cost, the idea is to consider a sensitivity matrix of lower 
dimensions. The second algorithm consists in combining the finest subdivi- 
sion with localization in time. The number of sections in this case coincides 
with one half of the number of parameters of the finest subdivision . At 

each iteration fc, for every section i = 1, . . . ,71^, Ug — — computing € 

W^y^^t costs n''g'^\ '01° where n^g'^^ denotes the cardinality 

of U O'*^. Moreover computing the SVD to obtain the new iteration has 

cost 4n2 (^^^-or^^ nf'''^ + s'^t^^ny(nf'^f + {nf'^ f. Finally computing 

the new prediction error has cost NNf^ . Although an higher number of systems 
must be solved, the algorithm is less costly since the sensitivity matrix has much 
lower dimensions. 
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Another possibility to decrease the cost of algorithm one, is to use the third 
algorithm, which consists in adopting an adaptive parametrization. At the A:— th 
iteration computing the sensitivity matrix G K"!'^^"» has cost rSg'^ NNf^, 
where the number of parameter n'g^ varies during the iterations and n'g^ < rig^ . 
Moreover computing the SVD to obtain the new iteration has cost AriyN'^nf'^ + 
9)Nny{n^g^Y + ^i^^e^Y- Finally computing the new prediction error has cost 
NN^. The gain with respect to the first strategy is evident if Ug'^ << rig^K 

The fourth algorithm combines both time localization and the adaptive 
parametrization. The number of sections in this case coincides with one half 
the number of parameters of the initial coarse subdivision 5^^^. The differ- 
ence with respect to the second algorithm is that the number of sections Us is 
lower, because it is no more related to the finest subdivision: in fact the adap- 
tive parametrization guides the choice of parameters to be estimated at each 
iteration. However the introduction of the adaptive parametrization introduces 
an inner loop. In detail, at each iteration fc, for every section i = 1, . . . 
applying the adaptive procedure until a minimum is reached (index I), com- 

puting ^^f'•^'''^'> e R"« ^ nt" xri'"-'-" ^^^^^ nf''''^\ ^' )Nl. Moreover com- 

puting the SVD to obtain the new iteration has cost ( ^ j,^" j rig + 

i^—jy^nyirSg'^'^^Y' -f [rSg''"'^^)'^. Finally computing the new prediction error 
has cost NNl- 

Just to give an idea of the computational gain of the fourth algorithm, 
computational costs of the four algorithms are summarized in table [2| averaging 
results of tests presented in section |6.4[ 





Finest subdivision 


Finest subdivision 
-|- time localization 


Adaptive subdivision 


Adaptive subdivision 
-|- time localization 


Computational 
cost 


8 ■ Ifll* 


2 . 10l3 


5 . 10l2 


S ■ 10^^ 



Table 2: Estimated computational cost of the four algorithms: using the finest subdi- 
vision and the projected damped Gauss Newton method, using the finest subdivision 
and the localization in time, using the adaptive parametrization and using the adaptive 
parametrization and time localization. 



6.4- Numerical results 

In this section we present some numerical tests to verify the effectiveness of 
the algorithm. As in section [5. 3 [ experimental data are simulated numerically, 
on 17 = [0, 8] X [0, l],Th = [0, 8] x {1} U [0, 8] x {0}. Moreover the velocity field 
u is modeled as a Poiseuille flow i.e. 

u(a:i,a;2)=( j. 
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We assume that ^ = 0.1, a = 0.1 and Cup ~ 0.1. Moreover we consider 
the finest subdivision with step length Aa; = 0.5. In algorithm 2 we con- 
sider {^1, . . . coincident with the finest subdivision, while in algorithm 4 
ns = 2 and {^i, . . . , = {0, 4, 8}. Define optimal subdivision the one which 
describes the real profile with the minimum number of parameters using the bi- 
section criterium. With distance from the optimal subdivision we indicate the 
number of points added (sign -|-) or subtracted (sign -) to the optimal sundivi- 
sion. We consider 9 test cases: results for the adaptive strategy with localization 
in time are shown in figure [9j In table [3] four algorithms are compared: us- 
ing the finest subdivision and the projected damped Gauss Newton method, 
using the finest subdivision and the localization in time, using the adaptive 
parametrization and using the adaptive parametrization and time localization. 

First of all observe that the number of iterations of algorithms 2 and 4 is 
higher since also sub-iterations to reach the minimum inside each section are 
counted (inner loop). 

In tests 1, 2 and 7, also working on the finest subdivision performs well, but 
it is much more costly. When the condition number of the sensitivity matrix 
increases, the accuracy is low. In particular in tests 3, 4, 5 it is evident 
how time localization improves convergence results both in algorithms 2 and 
4, selecting only some rows of However adopting only time localization 

is not sufficient in tests 6,7,9. Using an adaptive parametrization corresponds 
to select only some columns of ^'^j, considering a less number of parameters: 
in algorithm 3 the number of points added to the optimal subdivision is very 
low, but in general the estimates tend to be too much approximated. The 
best strategy consists in combining both adaptive parametrization and time 
localization (algorithm 4): this is a good compromise between good estimates 
and reasonable computational cost. Its effectiveness is evident e.g. in tests 8 
and 9. Moreover it only adds few points to the optimal subdivision. 



In figure 10 and 11 different iterations of the algorithm are shown for test 
8: it is evident how the algorithm firstly optimize parameters of section S2 = 
[4, 



X [0,1] (figure 11), and then that of si = [0,4] x [0,1] (in figure 10 the 
first 7 iteration are identical to the first one, thus only iterations 1 and 8 are 



plotted). The estimated subdivision is sketched in figure 12 it is evident how 
algorithm 4 slightly over-refine the optimal subdivision. 

6.5. Conditioning of the problem 

The ill-conditioning of the system matrix could increase when smaller 
segments are considered in Th : in fact in this case consecutive columns tend to be 
close to linear dependence, due to the small distance (Aa;) of the corresponding 
nodes in . This can be demonstrated numerically: consider in fact example 2 
presented in section |5.3| and generalize it considering the following parametric 
problem 



F„ ^[5-h,5]x {1} U [2 - /i, 2] X {0} , -d = (100, 80), <h<2. 
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Even supposing to know source localition Tin, solving the problem for different 
values oi h = {0.0625,0.125,0.25,0.5,1,2} and computing the condition num- 
ber of the sensitivity matrix, it can be seen that as h decreases, the condition 
number increases (cfr. figure 13). Since the condition number of the sensitivity 
matrix could become higher when smaller segments are considered, working on 
the finest subdivision could not be effective to reduce the ill-conditioning of the 
problem and an adaptive parametrization should be preferred. Observe more- 
over that in adaptive algorithms the Gauss Newton method is applied only to 
those parameters belonging to A'^'"'-': avoiding parameters less than the threshold 
62 is useful to reduce columns linear dependence. 

Moreover, as analyzed in section [6. 1[ at the stationary regime, the problem 
becomes ill-conditioned: thus, considering only the transitional regime, time 
localization could limit the ill-conditioning of the problem. This is evident e.g. 
in figure [Mj where the four algorithms are compared: without time localization 
(red dotted line) the condition number of the sensitivity matrix has a much 
higher upper bound. 



6. 6. Sensitivity of the fourth algorithm to thresholds variations 

It is interesting to analyze what happens when thresholds used in the fourth 
algorithm are changed, ei decides when a the segment corresponding to a pa- 
rameter should be refined: it is important to keep it not too low, to avoid 
over-refinements. £2 is such that parameters less than it are not considered to 
build the sensitivity matrix: avoiding small parameters reduces computational 
cost and the ill-conditioning of the problem, since we expect that they are not 
effective in output variations. 

Previous observations are summarized in table |4j where test 1 is consid- 
ered to understand how convergence results varies when thresholds are slightly 
changed: when ei is decreased the over-refinement increases, while when £2 is 
lower both the computational cost (number of iterations) and the condition 
number increase. When both ei and £2 decrease both the distance from the 
optimal subdivision and the computational cost and the condition number in- 
crease. Thus in general to reduce the cost is it better to increase £1, while to 
obtain more accurate results it could be useful to adopt smaller £1 and £2. 



7. Conclusions 

This paper presents a mathematical algorithm to solve a class of parabolic 
inverse problems based upon a convection-diffusion-reaction equation, extending 
some ideas presented in [8] and [22]. Both liquid (e.g. water) and gas (e.g. 
air) pollution problems could be considered: when source location is known, 
we have demonstrated that the problem is well-posed and can be solved e.g. 
using the Projected damped Gauss Newton method. When P™ is unknown, we 
have proved that adaptive parametrization with time localization is an effective 
strategy to estimate a sparse vector of parameters. 
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It could be interesting to introduce also an unrefinement strategy, trying 



to get closer to the optimal subdivision. For example consider figure 15 the 



optimal strategy would estimate only one parameter in [1, 2] x {1}, and it would 
not bisect the segment [1,2]. Instead algorithm 4 bisects [1,2]: the problem 
here is that the direction of the convective field u produces an overestimate of 
the right hand side parameter of [1,2] and an underestimate of the left hand 
side one. 

Another interesting aspect could be the generalizzation of the problem to 
time varying boundary conditions on and to analyze more deeply the prob- 
lem when space-time varying velocity fields are considered. 
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Appendix A. The importance of stabilizing the problem 

Dealing with convection dominated problems (||u|| >> /i) could be prob- 
lematic, due to spurious oscillations caused by the standard FE method. The 
simplest way to stabilize the problem is to refine the mesh, i.e. to consider a 
higher number of degrees of freedom; otherwise on a coarse mesh a stabilization 
method such as SUPG, DW or GLS, to mention only some of them, should be 
used. To simplify the problem in the following we apply the simplest strategy, 
i.e. we refine the mesh. However a stabilization method could be included in 
the model, modifying the weak FE formulation. Stabilization techniques are 
used e.g. in [H 

In this section we want to point out that the problem must be stabilized to 
obtain a correct estimate. In fact consider — [0, 8] x [0, 1], = [0, 8] x {1} U 
[0, 8] X {0}, the velocity field u is modelled as a Poiseuillc flow i.e. 

, . f -'^vxl + AVX2 \ 

u(a;i,a;2) = ( I , 

assume moreover that fi = 0.1, cr = 0.1, Cup — 0.1 and Fi„ = [0.5,1] x {1}, 
= 100. Apply to it the adaptive strategy with time localization, on different 
meshes. Results are depicted in figure |AT6) As it can be seen, when the mesh is 
too coarse, the presence of spurious oscillations compromise the convergence of 
the algorithm to the real profile, whereas adopting a fine mesh eliminates them 
and gives a good estimate of the boundary control. 
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Figure 9: Ntne test cases: results of the adaptive strategy with time localization: com- 
puted estimate (blu dotted line), real control (red line). For each figure: cost function 
(first row, left), error (first row, right), approximation of the upper horizontal seg- 
ment (second row, left), approximation of the bottom horizontal segment (second row, 
right). 
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Table 3: Comparison between four algorithms: using the finest subdivision and the pro- 
jected damped Gauss Newton method, using the finest subdivision and the localization 
in time, using the adaptive parametrization and using the adaptive parametrization 
and time localization. -error in the upper and lower horizontal segments, number 
of points added to the optimal subdivision in the upper and lower horizontal segments, 
number of iterations and final cost function. 
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Table 4: Test 1: results for different values of ei and €2. 
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Figure 10: Test 8. Adaptive parametrization and time localization. Evolution of the 
approximation (blu dotted line), real control (red line). Upper horizontal segment. 



Figure 11: Test 8. Adaptive parametrization and time localization. Evolution of the 
approximation (blu dotted line), real control (red line). Bottom horizontal segment. 
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Figure 12: Test 8. First row: optimal subdivision that could be obtained using a bisection 
strategy. Second row: estimated subdivision. 
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Figure 13: Example 2, with Ti„ = [5 — /i, 5] X {1} U [2 — /i, 2] X {0}. Condition number 
o/* ; for different values of h = {0.0625,0.125,0.25,0.5, 1,2}. 




Figure 14: Condition number of the sensitivity matrix with (blue line) and without (red 
dotted line) time localization. First row: finest subdivision. Second row: adaptive 
parametrization. Left: test 3. Right: test 9. 




Figure 15: Need of an under-refinement strategy. 
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Figure A. 16: Importance of using stabilization: concentration field (left), estimated 
profile (right). First row: using 41 nodes along x-axis and 9 along y-axis. Second row: 
using 81 nodes along x-axis and 13 along y-axis. Third row: using 81 nodes along 
X-axis and 21 along y-axis. 
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